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Abstract 

Monte Carlo simulation has been performed on a classical two dimensional XY- 
model with a modified form of interaction potential to investigate the role of topological 
defects on the phase transition exhibited by the model. In simulations in a restricted 
ensemble without defects, the system appears to remain ordered at all temperatures. 
Suppression of topological defects on the square plaquettes in the modified XY-model 
leads to complete elimination of the phase transition observed in this model. 

PACS: 05.10.Ln, 03.75.Lm, 64.60.an 



1 Introduction 

In 1984, Domany, Schick and Swendsen [1] introduced an extension of the two-dimensional 
(2D) XY-model where the classical spins (of unit length), located at the sites of a square 
lattice and free to rotate in a plane, say the XY plane (having no z-component) interact 
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with nearest-neighbors through a modified potential 
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where is the angle between the nearest neighbor spins and p 2 is a parameter used to 
change the shape of the potential. For p 2 = 1, the potential reduces to that of a conven- 
tional XY- model while with the increase in p 2 , the potential well gets narrower with a 
width ~ 7t/p and for 9 > n/p it is essentially constant at V(ir) = 2. The shape of the 
potential is shown in Fig. [T]for several values of p 2 . The conventional 2D XY-model does 
not possess any true long range order [2] and it is known that in this model, non-singular 
spin wave excitations alone cannot destroy the quasi long range order (QLRO). However, 
the presence of topological defects leads to a QLRO-disorder phase transition, now familiar 
as the Kosterlitz-Thouless (KT) phase transition. Kosterlitz and Thouless [31 H]predicted 
that topological singularities are both necessary and sufficient for the QLRO-disorder phase 
transition in the 2D XY-model and using a renormalization group (RG) approach they es- 
tablished that the phase transition is mediated by unbinding of vortices and anti-vortices 
which are stable topological defects in this system. The phase with QLRO is characterized 
by a slow algebraic decay of the spin-spin correlation function whereas a fast exponential 
decay is observed in a disordered system. The KT transition in the 2D XY model was 
unambiguously confirmed numerically by A. C. Irving and R. Kenna [5]. 

The modified XY-model of Eqn. (rjQ) has been analyzed by a number of investigators 
PQ El LTl El El HSl HU H2l H31 Ull US] and they all were of the opinion that it exhibits a 
first order phase transition for large values of p 2 . However, some investigators [161 El [18] 
attempted to interpret the Monte Carlo (MC) results for the first order phase transition 
in other ways. In a recent work [19J we have shown by performing extensive numerical 
simulations on relatively larger lattice sizes (upto 192 x 192) that the modified XY-model 
for large values of p 2 exhibits first order phase transition and all the finite size scaling 
rules for a first order phase transition were seen to be obeyed accurately, van Enter and 
Shlosman [121 [20] provided a rigorous proof of a first order phase transition in various SO(n)- 
invariant n-vector models which have a deep and narrow potential well and the model under 
investigation is a member of this general class of systems. 

It was argued in Ref. [T] that more than one type of excitations (i.e., topological exci- 
tations as well as vacancy excitations) may play a role in changing the nature of the phase 



transition, as one would become relatively more important than the other by the alteration 
of the potential. Later Himbergen argued that topological excitations alone are sufficient 
to account for both the continuous and the first order phase transitions, but in a qualita- 
tively different manner [5]. Algebraic topology or homotopy theory in the study of defects 
[2T| 122] has a wide application in the physics of phase transition. In the present paper, we 
investigate the role of topological defects in the phase transition of the modified XY-model 
under consideration. We are specifically interested in inquiring whether the first order phase 
transition in the modified XY-model is defect driven or not. In other words, if in the ab- 
sence of the role played by topological defects, would one observe the same order- disorder 
phase transition as the one found in the system with topological defects? If suppression of 
the defects changes the nature of the phase transition or eliminates it altogether, one may 
conclude that topological defects are necessary to describe the phase transition correctly. 

The present work was motivated to a great extent by the work of Lau and Dasgupta [23] 
and Dutta and Roy [24]. Lau and Dasgupta have shown that hedgehogs (point singulari- 
ties in the 3D Heisenberg model) are necessary for the phase transition in 3D Heisenberg 
model. They observed that if the formation of topological defects is suppressed in the 3D 
Heisenberg model, the system remains ordered at all temperatures and the transition to 
the disordered phase disappears altogether. Dutta and Roy performed MC simulations on 
planer P2 and P4 models which are known to exhibit a continuous and a first order phase 
transition respectively. It was noticed that the suppression of the defects in these models 
leads to a total disappearance of the phase transitions observed in these systems [23] . Other 
work, along the same line, that should be mentioned in this context, is that of Lammert 
et.al [25]. These authors have shown in a MC study that the nature of the nematic-isotropic 
transition changes when one suppresses the formation of the stable line defects, called the 
disclination lines. 

In the present work, we have found that topological defects are necessary for the existence 
and also for proper description of phase transitions exhibited by these class of models. We 
also find that the change in the nature of the phase transition that is observed with the 
change in the value of p 2 is due to a change in the role played by the topological defects in 
these systems. 

We arrange the rest of the paper as follows : Section 2 describes the simulational proce- 



dures used in the present work. The results and discussions are then presented in Section 3 
followed by the concluding remarks. 

2 The simulation details 

In order to study the behavior of the topological excitations and the role of topological 
defects in the phase transitions exhibited by the model under investigation, we have used 
the conventional Metropolis single spin update algorithm [26| [27] with some modifications 
in our MC simulation. We have found that while simulating a continuous lattice spin 
model using standard Metropolis algorithm, we need to adjust a parameter very carefully 
to generate a new a configuration. This parameter determines the amplitude of the random 
angular displacements of the spins and the results become very sensitive to the value of 
this parameter. In order to get rid of this difficulty of choosing the parameter, we generate 
a new spin configuration following the prescription of Wolff [28]. We take a random unit 
vector r and a spin flip s — > s' is defined as s ' = s — 2 (s, r)r where (s,r) is the dot 
product of s and r. Apart from this method of generating a new configuration, the rest of 
the algorithm is the standard Metropolis algorithm. Defining a spin flip in that way, our 
modified Metropolis algorithm is free from tuning any adjustable parameter while simulating 
a lattice spin model with continuous energy spectrum while the conditions of ergodicity and 
detail balance remain fulfilled. The Metropolis algorithm runs as follows: first we choose a 
spin at random from within the specified range. Then the change in energy AE associated 
with an attempted move is calculated. If AE < 0, the attempted move is accepted. If 
AE > 0, the attempted move is accepted with probability exp(— AE/T) where T is the 
dimensionless temperature. 

The average defect pair density is calculated in the following way. A vortex (anti-vortex) 
is a topological defect in which the angle variable 9, specifying the direction of the order 
parameter, changes by 2n(— 27r) in one circuit of any closed contour enclosing the excitation 
core. In order to trace out the topological defects, we consider a square plaquette in the 
physical space. Let s\, s*2, S3 and S4 be the four spins at the corners of the square plaquette. 
The angles between these adjacent spins are calculated with proper sign and these are then 
summed algebraically to find the total angle. The square plaquette is said to enclose a 



vortex (topological charge Q — 1) when the sum equals 2rr or more precisely very close 
to 2tc, taking into account the possible numerical errors. The square plaquette is said to 
enclose a anti-vortex (topological charge Q = — 1) if the sum equals — 2n. If the sum is zero, 
there is no topological defect in the plaquette. Average defect pair density (taking into 
consideration both vortices and anti-vortices) is calculated as the thermodynamic average 
of the absolute value of the vorticity summed over the entire lattice divided by the total 
number of spins. In this method, it is ensured that the net topological charge is always equal 
to zero in a system with periodic boundary conditions. It should be mentioned here that the 
smallest part of the system in real space that enclose a Q = ±1 point defect is a triangle. 
We can divide each elementary square plaquette diagonally into two triangles. One could 
thus consider a triangular plaquette in the physical space as well to trace out topological 
defects. We have tested that the total number of topological charges in the entire lattice 
remains the same whether we choose a square plaquette or a triangular plaquette. Only 
topological charges of strength Q = ±1 are considered since they are energetically favorable. 
In our investigation of the equilibrium behavior of topological defects near phase transition, 
we carried out simulations on system sizes with linear dimension L = 16, 32, 48 and 64 with 
periodic boundary conditions. 10 5 Monte Carlo sweeps (MCS) were used for equilibration 
and 10 6 MCS were used for calculating thermodynamic averages. One MC sweep is said to 
be completed when the number of attempted single spin moves equals the total number of 
spins in the system. The values of p 2 taken to study the variation of average defect pair 
density with p 2 are 4, 9, 16, 25, 36, 50, 64, 81 and 100. 

In order to implement the procedure of the suppression of topological defects in our 
model, a "chemical potential" term associated with the topological charges is included [231 
124"] . The modified Hamiltonian in the simulation is given by 



where 9ij is the angle between the nearest neighbor spins i, j and \Qijki\ is the absolute 
value of the charge enclosed by the square plaquette. A positive value of A ensures that 
the formation of the charges becomes expensive in terms of energy. So for positive A, this 
"chemical potential" term has the effect of suppressing configurations containing defects. In 
simulation with this modified Hamiltonian, a calculation of AE, the energy change associ- 





ated with an attempted move of a spin, involves calculations of the changes in topological 
charges associated with the four unit square plaquettes which share the spin under consid- 
eration. For almost complete suppression of the defects, the value of A was chosen to be 
between 5 and 20 irrespective of temperature. The A — > oo limit of Eqn. ([2]) indicates an 
ensemble in which configurations containing topological defects are not allowed. We started 
our simulation in that restricted ensemble with a configuration in which all the spins are 
aligned parallel to one another, ie, there is no topological defect. The restricted simulations 
were carried out by using the modified Metropolis spin update algorithm described earlier in 
this section. We performed our restricted simulations on system sizes with linear dimension 
L = 16, 32, 64 and 96 with periodic boundary conditions. The A = corresponds to an 
unrestricted simulation, where no suppression of topological defects take place. 

For the purpose of calculating various thermodynamic quantities, we have used multiple 
histogram reweighting technique of Ferrenberg and Swendsen [29]. In the restricted simula- 
tions, 10 6 MCS were taken for equilibration and 10 7 MCS were used for computing the raw 
histograms (both energy and order parameter histograms). The value of p 2 was taken to be 
50 in order to carry out the restricted simulations. 

3 Results and discussions 

In this section we present in detail the results obtained from our simulations. 

3.1 Behavior of topological excitations near phase transition 

We used the method described in Sec. [2] to determine the average defect pair density (p) 
of the system. The variation of the p with the dimensionless temperature T is shown in 
Fig. [2] for several values of the parameter p 2 . The average defect pair density is found 
to increase sharply as T increases through the transition temperature T c (p 2 ) and appears 
to exhibit a sharp jump at T c (p 2 ), particularly for p 2 ^> 1. Fig. [2] indicates that T c (p 2 ) 
decreases as the values of p 2 increases. It is evident from Fig. [2] that for larger values of p 2 , 
at some temperature T c (p 2 ), vortices suddenly appear in great numbers and a first order 
phase transition takes place which is in accordance with the explanation of Himbergen [6]. 
It may be noted that Jonsson, Minnhagen and Nylen [11] also performed MC simulations 



on a 2D XY-model with a modified potential, which is essentially equivalent to that of Eqn. 
dlj), and interpreted the first order transition to be of vortex unbinding type. 

We have also studied the behavior of topological excitations with the parameter p 2 . The 
average defect pair density (p) as a function of the parameter p 2 is plotted in Fig. [3] for three 
different system sizes at a temperature T = 1.12 which is above the transition temperature 
of the model for p 2 = 50. We observe that above the transition temperature, the data for p 
versus p 2 are nicely fitted by the following expression 

P(T) = Pmax - a(T)exp(-7Vp 2 ) (3) 

Eqn. ([3]) takes into account both vortices and anti-vortices. The values of p max for the three 
system sizes are listed in Table 1. There is no significant system size dependence of the 
parameters and it may be noted that p increases with p 2 . In the limit p 2 — > oo, the system 
contains only vortex excitations. This means that in the high p 2 limit, the system must 
be disordered even at very small temperatures and consequently the transition temperature 
must be very low. This is the reason behind the decrease in the transition temperature with 
increase in p 2 . In this context, we refer to the work of Romano et.al [30J who, in a MC 
study of 2D generalized XY-model with spin component n = 3, discussed the nature of the 
phase transition with the variation of a generalized parameter. 

We have also calculated the defect core energy (E c ) for various values of p 2 . Because 
there is always a positive energy cost E c associated with the creation of a vortex core, 
thermally excited vortices in thermal equilibrium always contribute terms proportional to 
exp(—E c /T) to the partition function. Therefore the total number of topological charges 
(n) shows an exponential behavior exp(— E c /T) at low temperatures and is given by 

n = n exp(-E c /T) (4) 

Taking natural logarithm on both sides of Eqn. (j3J) 

E c 

In n = In n - — (5) 

The defect core energy E c for each p 2 is determined from the linear fit of the plot In n versus 
1/T and Fig. H] shows the plot of variation of E c versus p 2 . The defect core energy E c for 
p 2 = 50 model (which is known to exhibit a strongly first order phase transition) is found 



to be 11.911 ± 0.15 while that for p 2 = 1 model (which is known to exhibit a continuous 
phase transition) is found to be 7.560 ± 0.019. This is in apparent contradiction with the 
main finding of Saito [31] who, in a MC study of a system of interacting dislocation vectors, 
predicted a continuous phase transition due to large core energy and a first order phase 
transition due to small core energy. 

3.2 Restricted simulations with no defects 

Before presenting our results of restricted simulations, we briefly define the thermodynamic 
quantities that we have evaluated. The MC simulations were carried out with the modified 
Hamiltonian given by Eqn. (J2]) where the new term acts as a "chemical potential" for the 
defects. 

The specific heat C v is evaluated from the energy fluctuations 

c ( 6) 

where T is the dimensionless temperature and iV = L 2 is the total number of spins. The 
average order parameter is given by 

(Pi) = (cos0) (7) 

where <fi is the angle that a spin makes with the preferred direction of orientation and the 
average is over the entire sample. The first rank pair correlation function is defined as 

G 1 (r) = <(co B e ij )> r (8) 

where % and j are two spins separated by a distance r. The second rank pair correlation 
function is defined as 

G 2 (r) = (P 2 (cosd ij )) r (9) 

We did not find any evidence for a phase transition from the ordered to the disordered 
phase at any temperature in the simulations of the restricted ensemble where configurations 
containing defects are not allowed. 

The energy histograms obtained for L = 64 are shown in Fig. |5j For this lattice, 
simulations were performed at 12 different temperatures ranging from T = 1.0375 to T = 
1.1600. It is evident from the energy histograms that the dual peak nature of the histograms 



as obtained for an unrestricted simulation [TH] disappears in a restricted simulation. It is 
known that the dual peak nature of the histograms is a signature of a first order transition 
where two phases can coexist at a given temperature. 

Fig. O shows the temperature dependence of the average energy (E) for a number of 
lattices, as obtained by applying the histogram reweighting technique. For comparison, the 
same plots for unrestricted simulations (where no defect is suppressed) are shown by thick 
lines in the same figure. It is evident from the figure that the energy changes only gradually 
and smoothly with temperatures for the restricted simulations while a sharp variation of 
the same with temperature is observed in the unrestricted case. The average value of the 
order parameter (Pi) defined in Eqn. ([71) , is, always nonzero for a finite size system. Hence 
we have studied the system size dependence of (-Pi). The values of (Pi) at T — > oo versus 
1/L are plotted in Fig. Hand the system size dependence is, in fact, well fitted by the form 
(Pj) = P + a/L with P = 0.587 ± 0.005 and a = 3.035 ± 0.141. It is clear from Fig. [7J that 
there is no indication of (Pi) extrapolating to zero in the thermodynamic limit L — > oo and 
thereby suggesting a state with long range ferromagnetic order. 

The specific heat C v was obtained from the energy fluctuation relation (Eqn. ([6])). The 
specific heat data for L = 64 in a restricted ensemble are shown by dashed line in Fig.[8]where 
the results for the unrestricted case are also shown by solid line for comparison. For clarity 
the data for the restricted and the unrestricted simulations are plotted in two different scales. 
While C v has a large peak height (~ 700) at the transition temperature in the unrestricted 
case, which presumably is a signal of a phase transition in a finite system, in the defect-free 
case the peak height (~ 16) is drastically reduced and almost disappears in comparison 
with the normal case (where no defect is suppressed). We would like to argue that, in 
the restricted ensemble, the existence of a peak in C v of insignificant height (compared to 
that of an unrestricted ensemble) over the temperature range cannot be a sign of a phase 
transition. These may be attributed to the fact that complete suppression of topological 
defects is never possible, there always exists a small number of residual charges in the system. 
We have also calculated the free energy like quantity A from the energy histograms. It is 
defined as A(E; (3, L,J\f) = — \nN(E; (3, L) where N(E; (3, L) is the histogram count of the 
energy distribution. Fig. [9] shows the plot of the quantity A against E for L = 64. The 
inset of Fig. [9] shows the same plot for the original model (Eqn. ([I])) where a double well 



structure of equal depth at the transition temperature signals a first order transition. We 
observe the absence of any such double well structure in A when defects are suppressed. 
We would be inclined to conclude from the results of C v and A that the defect free phase 
exhibits no phase transition at all. 

We now turn to pair correlation functions defined earlier in this section. Fig. [10] shows 
the plot of Gi(r) against r for L = 64 in the restricted as well as the unrestricted cases. 
The first rank pair correlation function G\(r) for p 2 = 50 at temperatures T = 1.0500 and 
T = 1.1000, which is much higher than the transition temperature of the original model, 
decays exponentially to zero in the unrestricted simulations, as it should, in the complete 
absence of long range order and a best fit with G\(r) = aexp(— 5r) yields a = 0.745 ± 0.01 
and 5 = 0.700 ± 0.01 for T = 1.0500 and a = 0.806 ± 0.009 and 5 = 0.944 ± 0.009 for 
T = 1.1000. For the simulation where the defects are suppressed, Gi(r) decays algebraically 
and a best fit with G\(r) = ar~ b +f yields the parameter a = 0.064±0.003, b = 0.175±0.016 
and / = 0.806 ± 0.001 for T = 1.0500 and a = 0.019 ± 0.009, b = 0.271 ± 0.018 and 
/ = 0.590 ± 0.001 for T = 1.1000. It may be noted that the parameter / is the asymptotic 
value of the pair correlation function. The next higher order correlation function G^O") 
against r for L = 64 is plotted in Fig. [11] at T = 1.0500 and T = 1.0100 for both the cases. 
The results indicate that G^O") decays algebraically in both the cases and long range order 
prevails in the system via higher order correlation functions. 

We now need to address the question of phase space connectivity before arriving at the 
conclusion that topological defects are indeed necessary for the phase transition. Since we 
have used large values of A in our restricted simulations in order to suppress the evolution 
of topological defects, we have to demonstrate that the observed behavior is not caused by 
trapping of the system in a small region of phase space with nonzero (-Pi). Any MC study 
is guaranteed to generate appropriate ensemble averages if there is a path connecting any 
two points in the phase space with nonzero transition probability. We have investigated the 
phase space connectivity by observing the evolution of the order parameter and energy with 
MC sweeps. The connectedness is satisfied if the observed quantities for different initial 
states converge to the same final value. In Fig. [T2] we have shown that for L = 64, after 
suppressing the defects (by using A = 20) on the square plaquettes, the final values of the 
order parameter is same for three different initial configurations. This observation ensures 



that we can use a value of A upto 20 without violating the phase space connectivity and the 
observed non- vanishing of (Pi) is not a result of trapping of the system in the phase space. 

4 Conclusions 

It is established in this paper that topological defects play a very crucial role in the phase 
transitions exhibited by the models we discussed. We have observed that the average defect 
pair density grows rapidly with the increase in p 2 (which increases the non-linearity of 
the potential well). For high p 2 , the potential well becomes narrower so that there is 
an insufficient increase in the defect density at low temperatures and then, at a certain 
temperature, they suddenly appear in the system in great numbers. Therefore it may be 
thought that for larger values of p 2 the class of models, we have investigated, behave like 
a dense defect system and give rise to first order phase transition as has been predicted by 
Minnhagen [SJ |9j [10]. It has also been observed that the first order transition is eliminated 
totally when configurations containing topological defects are not allowed to occur and 
the system appears to remain ordered at all temperatures. Hence topological defects are 
necessary to account for the first order phase transition for larger values of p 2 . 

Another point which must be mentioned before ending this section is the performance 
of the Metropolis algorithm with the modification discussed earlier in this paper. The mod- 
ification makes us free from tuning any adjustable parameter while simulating a continuous 
model and this has resulted in the model being successfully simulated. 
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Table 1: parameters for the fit of p(T) = p max — a(T)exp(— 7^/p 2 ) for different L 
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Pmax 


a(T) 
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32 


0.2866 ± 0.002 


0.411 ±0.005 


0.301 ±0.008 


48 


0.2876 ± 0.002 


0.406 ± 0.006 


0.297 ±0.009 


64 


0.2878 ± 0.002 


0.406 ± 0.006 


0.296 ±0.009 
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Figure 1: The potential function of Eqn. (1) is shown for different values of p 2 . 




gure 2: Average defect pair density p plotted against dimensionless temperature T for 
= 64 for various values of p 2 
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Figure 3: Average defect pair density p plotted as a function of p 2 at T = 1.28 for three 
system sizes. The best fit corresponds to L — 64. The errorbars are smaller than the 
dimension of the symbols used for plotting. 
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Figure 4: The defect core energy E c plotted against the parameter p 2 with the linear fit 
represented by the dotted line. The error bars are of the dimension of the symbols used for 
plotting. 
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Figure 5: The histograms for E, the average energy per particle generated for the 64 x 64 
lattice for the p 2 = 50 model at the 12 temperatures indicated. 




Figure 6: The average energy per particle E plotted against dimensionless temperature T for 
different lattice sizes. The thick curves (on the left) correspond to unrestricted simulations. 
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Figure 7: The plot of order parameter (Pi) at T — 1.12 vs 1/L. The straight line is the best 
linear fit to the data. The errorbars are smaller than the dimension of the symbols used for 
plotting. 
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Figure 8: The specific heat C v plotted against dimensionless temperature T for restricted 
ensemble (represented by dashed line) for L = 64. The same is plotted for unrestricted 
ensemble (shown by solid line). For clarity the graphs are plotted in two different scales. 
The bottom X-axis and left Y-axis are chosen to plot the results for unrestricted simulations 
while the top X-axis and right Y-axis are chosen to plot the same for restricted simulations. 
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Figure 9: The free energy A generated by the multiple histogram reweighting technique 
plotted against energy per particle for L = 64 in the restricted ensemble. Absence of a 
double well structure is to be noted. The inset shows the same plot for an unrestricted 
ensemble where the presence of a double well structure is observed. 




I 


i i 


A 




T 
♦\ 

\\ 
:\ 


T=1.0500^» 

T=1 .1000 — > a 

T=1 .0500 (unrestricted) -> t 

T=1.1000 (unrestricted) -> ♦ 


A 

i 


-« » » »- ♦ ♦ ♦— » ♦ ♦ » » ♦ — 
I I 



30 



r 

Figure 10: The plots of the pair correlation function G 1 (r) against r for the 64 x 64 lattice 
for the temperatures indicated. The curves are plotted for r ranging upto L/2. 
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Figure 11: The plots of the pair correlation function G 2 (r) against r for the L = 64 lattice 
for the temperatures indicated. The curves are plotted for r ranging upto L/2. 




Figure 12: The evolution of the order parameter at T = 1.1800 for L = 64 lattice after 
suppressing the defects using A = 20 for three different initial configurations: Pi = 0.999, 
Pi = 0.244 and Pi = 0.012. The final values of the order parameter have the same average 
value. 



